﻿///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using System;
using ILNumerics;
using ILNumerics.BuiltInFunctions;

// Exponential functions
// =====================
//
// exp - Exponential 
// expm1-  Compute exp(x)-1 accurately for small values of x 
// log - Natural logarithm 
// log10 - Common (base 10) logarithm 
// log1p - Compute log(1+x) accurately for small values of x 
// log2 - Base 2 logarithm and dissect floating-point numbers into exponent and mantissa 
// nextpow2 - Next higher power of 2 
// nthroot - Real nth root of real numbers 
// pow2 - Base 2 power and scale floating-point numbers 
// power - Array power 
// reallog - Natural logarithm for nonnegative real arrays 
// realpow - Array power for real-only output 
// realsqrt - Square root for nonnegative real arrays 
// sqrt - Square root

namespace MathLib
{
    /// <summary>
    /// Exponential functions not provided by ILNumerics.
    /// </summary>
    public partial class MLMath
    {
        //
        // --- Compute nth root ---
        //

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<double> arg, double n)
        {
            return ILMath.pow(arg, 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<float> nthroot(ILArray<float> arg, double n)
        {
            return ILMath.pow(arg, (float)(1.0 / n));
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<int> arg, double n)
        {
            return ILMath.pow(ILMath.todouble(arg), 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<complex> nthroot(ILArray<complex> arg, complex n)
        {
            return ILMath.pow(arg, 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<fcomplex> nthroot(ILArray<fcomplex> arg, complex n)
        {
            return ILMath.pow(arg, (fcomplex)(1.0 / n));
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<byte> arg, double n)
        {
            return ILMath.pow(ILMath.todouble(arg), 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<char> arg, double n)
        {
            return ILMath.pow(ILMath.todouble(arg), 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<long> arg, double n)
        {
            return ILMath.pow(ILMath.todouble(arg), 1.0 / n);
        }

        /// <summary>
        /// Compute nth root.
        /// </summary>
        /// <param name="arg">Argument to return nth root of</param>
        /// <param name="n">Number of root to return</param>
        /// <returns>nth root of argument</returns>
        public static ILArray<double> nthroot(ILArray<short> arg, double n)
        {
            return ILMath.pow(ILMath.todouble(arg), 1.0 / n);
        }


        //
        // --- Compute exp(x) - 1 accurately for small x ---
        //

        /// <summary>
        /// Computes exp(x) - 1 accurately for small values of x.
        /// </summary>
        /// <param name="power">Power of e</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> expm1(ILArray<double> power)
        {
            ILArray<double> expm1Array = new ILArray<double>(power.Dimensions.ToIntArray());

            for (int i = 0; i < expm1Array.InternalArray4Experts.Length; i++)
            {
                double x = power.InternalArray4Experts[i];

                // exp(x) - 1 = x/2 + x2/6 + x3/24 + x4/120
                const double c1 = 1.0 / 2.0;
                const double c2 = 1.0 / 6.0;
                const double c3 = 1.0 / 24.0;
                const double c4 = 1.0 / 120.0;
                double x2 = x * x;
                double x3 = x2 * x;
                double x4 = x3 * x;
                expm1Array.InternalArray4Experts[i] = c1 * x + c2 * x2 + c3 * x3 + c4 * x4;
            }

            return expm1Array;
        }

        /// <summary>
        /// Computes exp(x) - 1 accurately for small values of x.
        /// </summary>
        /// <param name="power">Power of e</param>
        /// <returns>Complex array containing result</returns>
        public static ILArray<complex> expm1(ILArray<complex> power)
        {
            ILArray<complex> expm1Array = new ILArray<complex>(power.Dimensions.ToIntArray());

            for (int i = 0; i < expm1Array.InternalArray4Experts.Length; i++)
            {
                complex x = power.InternalArray4Experts[i];

                // exp(x) - 1 = x/2 + x2/6 + x3/24 + x4/120
                const double c1 = 1.0 / 2.0;
                const double c2 = 1.0 / 6.0;
                const double c3 = 1.0 / 24.0;
                const double c4 = 1.0 / 120.0;
                complex x2 = x * x;
                complex x3 = x2 * x;
                complex x4 = x3 * x;
                expm1Array.InternalArray4Experts[i] = c1 * x + c2 * x2 + c3 * x3 + c4 * x4;
            }

            return expm1Array;
        }


        //
        // --- Compute log(1 + x) accurately for small x ---
        //

        /// <summary>
        /// Computes log(1 + x) accurately for small values of x.
        /// </summary>
        /// <param name="arg">Value of x to compute log(1 + x) for</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> log1p(ILArray<double> arg)
        {
            ILArray<double> log1pArray = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < log1pArray.InternalArray4Experts.Length; i++)
            {
                double x = arg.InternalArray4Experts[i];

                // log(x + 1) = x - x2/2 + x3/3 - x4/4
                const double c2 = -1.0 / 2.0;
                const double c3 = 1.0 / 3.0;
                const double c4 = -1.0 / 4.0;
                double x2 = x * x;
                double x3 = x2 * x;
                double x4 = x3 * x;
                log1pArray.InternalArray4Experts[i] = x + c2 * x2 + c3 * x3 + c4 * x4;
            }

            return log1pArray;
        }

        /// <summary>
        /// Computes log(1 + x) accurately for small values of x.
        /// </summary>
        /// <param name="arg">Value of x to compute log(1 + x) for</param>
        /// <returns>Complex array containing result</returns>
        public static ILArray<complex> log1p(ILArray<complex> arg)
        {
            ILArray<complex> log1pArray = new ILArray<complex>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < log1pArray.InternalArray4Experts.Length; i++)
            {
                complex x = arg.InternalArray4Experts[i];

                // log(x + 1) = x - x2/2 + x3/3 - x4/4
                const double c2 = -1.0 / 2.0;
                const double c3 = 1.0 / 3.0;
                const double c4 = -1.0 / 4.0;
                complex x2 = x * x;
                complex x3 = x2 * x;
                complex x4 = x3 * x;
                log1pArray.InternalArray4Experts[i] = x + c2 * x2 + c3 * x3 + c4 * x4;
            }

            return log1pArray;
        }


        //
        // --- Compute logarithm base 2 ---
        //

        /// <summary>
        /// Computes logarithm base 2.
        /// </summary>
        /// <param name="arg">Value of x to compute log2(x) for</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> log2(ILArray<double> arg)
        {
            ILArray<double> log2Array = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < log2Array.InternalArray4Experts.Length; i++)
            {
                double x = arg.InternalArray4Experts[i];
                log2Array.InternalArray4Experts[i] = Math.Log(x, 2);
            }

            return log2Array;
        }

        /// <summary>
        /// Computes logarithm base 2.
        /// </summary>
        /// <param name="arg">Value of x to compute log2(x) for</param>
        /// <returns>Complex array containing result</returns>
        public static ILArray<complex> log2(ILArray<complex> arg)
        {
            ILArray<complex> log2Array = new ILArray<complex>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < log2Array.InternalArray4Experts.Length; i++)
            {
                complex x = arg.InternalArray4Experts[i];
                log2Array.InternalArray4Experts[i] = complex.Log(x) / complex.Log(2);
            }

            return log2Array;
        }


        //
        // --- Compute next power of 2 ---
        //

        /// <summary>
        /// Computes next power of 2.
        /// </summary>
        /// <param name="arg">Value of x to calculate next power of two for</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> nextpow2(ILArray<double> arg)
        {
            ILArray<double> np2Array = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < np2Array.InternalArray4Experts.Length; i++)
            {
                int x = (int)arg.InternalArray4Experts[i];
                int np2 = 0;
                while (x > 1)
                {
                    x >>= 1;
                    np2++;
                }
                np2Array.InternalArray4Experts[i] = np2;
            }

            return np2Array;
        }


        //
        // --- Compute power of 2 ---
        //

        /// <summary>
        /// Computes power of 2.
        /// </summary>
        /// <param name="arg">Value of x to compute pow2(x) for</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> pow2(ILArray<double> arg)
        {
            ILArray<double> pow2Array = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < pow2Array.InternalArray4Experts.Length; i++)
            {
                double x = arg.InternalArray4Experts[i];
                pow2Array.InternalArray4Experts[i] = Math.Pow(2, x);
            }

            return pow2Array;
        }

        /// <summary>
        /// Computes power of 2.
        /// </summary>
        /// <param name="arg">Value of x to compute pow2(x) for</param>
        /// <returns>Complex array containing result</returns>
        public static ILArray<complex> pow2(ILArray<complex> arg)
        {
            ILArray<complex> pow2Array = new ILArray<complex>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < pow2Array.InternalArray4Experts.Length; i++)
            {
                complex x = arg.InternalArray4Experts[i];
                pow2Array.InternalArray4Experts[i] = complex.Pow(2, x);
            }

            return pow2Array;
        }

        /// <summary>
        /// Computes power of 2.
        /// </summary>
        /// <param name="scale">Scaling factor</param>
        /// <param name="exp">Exponent</param>
        /// <returns>Real array containing result</returns>
        public static ILArray<double> pow2(ILArray<double> scale, ILArray<int> exp)
        {
            ILArray<double> pow2Array = new ILArray<double>(scale.Dimensions.ToIntArray());

            for (int i = 0; i < pow2Array.InternalArray4Experts.Length; i++)
            {
                double f = scale.InternalArray4Experts[i];
                int e = exp.InternalArray4Experts[i];
                pow2Array.InternalArray4Experts[i] = Math.Pow(2, e) * f;
            }

            return pow2Array;
        }

        /// <summary>
        /// Computes power of 2.
        /// </summary>
        /// <param name="scale">Scaling factor</param>
        /// <param name="exp">Exponent</param>
        /// <returns>Complex array containing result</returns>
        public static ILArray<complex> pow2(ILArray<complex> scale, ILArray<int> exp)
        {
            ILArray<complex> pow2Array = new ILArray<complex>(scale.Dimensions.ToIntArray());

            for (int i = 0; i < pow2Array.InternalArray4Experts.Length; i++)
            {
                complex f = scale.InternalArray4Experts[i];
                int e = exp.InternalArray4Experts[i];
                pow2Array.InternalArray4Experts[i] = Math.Pow(2, e) * f;
            }

            return pow2Array;
        }


        //
        // --- Raise matrix m elements to power of matrix p elements ---
        //

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<double> power(ILArray<double> m, ILArray<double> p)
        {
            ILArray<double> powerArray = new ILArray<double>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<float> power(ILArray<float> m, ILArray<float> p)
        {
            ILArray<float> powerArray = new ILArray<float>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (float)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (float)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<complex> power(ILArray<complex> m, ILArray<complex> p)
        {
            ILArray<complex> powerArray = new ILArray<complex>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                complex pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    complex me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = complex.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    complex me = m.InternalArray4Experts[i];
                    complex pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = complex.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<fcomplex> power(ILArray<fcomplex> m, ILArray<fcomplex> p)
        {
            ILArray<fcomplex> powerArray = new ILArray<fcomplex>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                fcomplex pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    fcomplex me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = fcomplex.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    fcomplex me = m.InternalArray4Experts[i];
                    fcomplex pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = fcomplex.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<int> power(ILArray<int> m, ILArray<int> p)
        {
            ILArray<int> powerArray = new ILArray<int>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<uint> power(ILArray<uint> m, ILArray<uint> p)
        {
            ILArray<uint> powerArray = new ILArray<uint>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (uint)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (uint)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<int> power(ILArray<byte> m, ILArray<byte> p)
        {
            ILArray<int> powerArray = new ILArray<int>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<int> power(ILArray<char> m, ILArray<char> p)
        {
            ILArray<int> powerArray = new ILArray<int>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (int)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<long> power(ILArray<long> m, ILArray<long> p)
        {
            ILArray<long> powerArray = new ILArray<long>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (long)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (long)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<ulong> power(ILArray<ulong> m, ILArray<ulong> p)
        {
            ILArray<ulong> powerArray = new ILArray<ulong>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (ulong)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (ulong)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<short> power(ILArray<short> m, ILArray<short> p)
        {
            ILArray<short> powerArray = new ILArray<short>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (short)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (short)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        /// <summary>
        /// Raises elements of matrix m to power of elements of matrix p.
        /// </summary>
        /// <param name="m">Base matrix</param>
        /// <param name="p">Power matrix</param>
        /// <returns>Array containing result</returns>
        public static ILArray<ushort> power(ILArray<ushort> m, ILArray<ushort> p)
        {
            ILArray<ushort> powerArray = new ILArray<ushort>(m.Dimensions.ToIntArray());

            if (p.InternalArray4Experts.Length == 1)
            {
                double pe = p.InternalArray4Experts[0];
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (ushort)Math.Pow(me, pe);
                }
            }
            else
            {
                for (int i = 0; i < powerArray.InternalArray4Experts.Length; i++)
                {
                    double me = m.InternalArray4Experts[i];
                    double pe = p.InternalArray4Experts[i];
                    powerArray.InternalArray4Experts[i] = (ushort)Math.Pow(me, pe);
                }
            }

            return powerArray;
        }

        //
        // --- Sqrt (Compute square root) ---
        //

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(double arg)
        {
            return Math.Sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<double> arg)
        {
            return ILMath.sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static float sqrt(float arg)
        {
            return (float)Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<float> sqrt(ILArray<float> arg)
        {
            return ILMath.sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(int arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<int> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(uint arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<uint> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static complex sqrt(complex arg)
        {
            return complex.Sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<complex> sqrt(ILArray<complex> arg)
        {
            return ILMath.sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static fcomplex sqrt(fcomplex arg)
        {
            return fcomplex.Sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<fcomplex> sqrt(ILArray<fcomplex> arg)
        {
            return ILMath.sqrt(arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(byte arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<byte> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(char arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<char> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(long arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<long> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(ulong arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<ulong> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(short arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<short> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static double sqrt(ushort arg)
        {
            return Math.Sqrt((double)arg);
        }

        /// <summary>
        /// Compute square root.
        /// </summary>
        /// <param name="arg">Argument to return square root of</param>
        /// <returns>square root of argument</returns>
        public static ILArray<double> sqrt(ILArray<ushort> arg)
        {
            return ILMath.sqrt(ILMath.todouble(arg));
        }
    }
}
